#PBAR_Zspec_CompareOldNewData.py
#
#  Compare old and new zspec data.
# 8/19/2013, John Kwong

import os
import numpy as np
import PBAR_Zspec
import matplotlib.pyplot as plt

# Location of data
basepath = r'C:\Users\jkwong\Documents\Work\PBAR\data'
acquisitionTime = [300., 300., 300., 300.]
pulseRate = [60., 60., 60., 60.]

countRange = np.array([0.7e3, 1e3])
countRange = np.array([0.7, 1])

#countRange = np.array([250, 325])

# range in which the spectra mush intersect with the horizontal range
binBounds = np.array([20, 200])

# detector to match
detStandard = 69

# list of good/bad detectors
(goodDetectorsList, badDetectorsList) = \
                    PBAR_Zspec.GenerateDefaultDetectorList()
# Load zspec data
filenameList = ['dx95.csv', 'ea64.csv', 'ea65.csv', 'ea68.csv']
filenameList = ['dx95.csv','ea68.csv']

fullfilenameList = []
for i in xrange(len(filenameList)):
    fullfilenameList.append(os.path.join(basepath, filenameList[i]))

dat = PBAR_Zspec.ReadZspec(fullfilenameList)

# get the bin values for all the datasets

# Shift among PMTs
binMeanCC = np.zeros((len(dat), dat[0].shape[1]))
gainShift = np.zeros((len(dat), dat[0].shape[1]))
#for i in xrange(len(dat)):
#    #bins = np.arange(dat[0].shape[0])
#    binArray = np.matlib.repmat(np.arange(0,dat[0].shape[0]),dat[0].shape[1],1).T
#    cutt = (binArray > binBounds[0]) & (binArray < binBounds[1]) & (dat[ii] > countRange[0]) & (dat[ii] < countRange[1])
#    binMean = (binArray * dat[ii] * cutt).sum(axis = 0).astype(float) / (dat[ii] * cutt).sum(axis = 0).astype(float) # min bin, calculated using the count as weight
#    binMean = binMean.astype(float)
#    pmtShift = binMean[detStandard] / binMean
#    pmtShift[np.isnan(pmtShift)] = -1
binArray = np.matlib.repmat(np.arange(0,dat[0].shape[0]),dat[0].shape[1],1).T.astype(float)
for ii in range(len(dat)):
    cutt = (binArray > binBounds[0]) & (binArray < binBounds[1]) & (dat[ii]/acquisitionTime[ii] > countRange[0]) & (dat[ii]/acquisitionTime[ii] < countRange[1])
    binMeanCC[ii,:] = (binArray * dat[ii] * cutt).sum(axis = 0).astype(float) / (dat[ii] * cutt).sum(axis = 0).astype(float)
    gainShift[ii,:] = binMeanCC[0,:] / binMeanCC[ii,:]


# Plot the waveforms at multiple times bins, SUBPLOTS

# index = 1
#index = 6
peakIndex = 0

plotFits = 1
filterSpectra = 0
filterWidth = 1
normalize = 0
correctGain = 0
plotRate = 1

subplotx = 4
subploty = 2

# cycle through detectors
figNum = 0
for detNum in np.arange(24,64):
    # Generate a new figure window if filled.
    if ((figNum % (subplotx * subploty)) == 0):
        f, ax = plt.subplots(subploty, subplotx, sharex='col', sharey='row')
    # Calculate subplot window index
    ii = (figNum % (subplotx*subploty)) /subplotx
    jj = (figNum % (subplotx*subploty)) % subplotx

    
    for i in xrange(1, 2):
        if correctGain:
            x = np.arange(256).astype(float) * gainShift[i,detNum]
        else:
            x = np.arange(256).astype(float)
        if plotRate:
            ax[ii][jj].plot(x, dat[i][:,detNum]/acquisitionTime[i] + 1e-10, label = '%s' %filenameList[i])
            ax[ii][jj].plot(x, dat[i][:,detNum]/acquisitionTime[i]*2.0 + 1e-10, label = '%s' %filenameList[i])
            ax[ii][jj].plot(x, dat[i][:,detNum]/acquisitionTime[i]*0.5 + 1e-10, label = '%s' %filenameList[i])
        else:
            ax[ii][jj].plot(x, dat[i][:,detNum] + 1e-10, label = '%s' %filenameList[i])

#        ax[ii][jj].plot(x, dat[i][:,detNum]/acquisitionTime[i]/pulseRate[i] + 1e-10, label = '%s' %filenameList[i])
    
    ax[ii][jj].grid()
    if (ii == 0) & (jj == 0):
        ax[ii][jj].legend()
    #ax[ii][jj].set_xlim((100, 200))
    ax[ii][jj].set_xlim((60, 160))
    
    
    ax[ii][jj].set_yscale('log')
    if plotRate:
        ax[ii][jj].set_ylabel('Rate')
        ax[ii][jj].set_ylim((1e-2, 1e1))
    else:
        ax[ii][jj].set_ylabel('Count')
        ax[ii][jj].set_ylim((1e-2, 1e6))
    ax[ii][jj].set_xlabel('Bin')
    ax[ii][jj].set_title('Det# %d' %(detNum+1))
    figNum += 1
    
# Plot ratio of waveform amlitudes

# index = 1
#index = 6
peakIndex = 0

plotFits = 1
filterSpectra = 0
filterWidth = 1
normalize = 0
correctGain = 0


subplotx = 4
subploty = 2

# cycle through detectors
for detNum in np.arange(80,88):
    # Generate a new figure window if filled.
    if ((detNum % (subplotx * subploty)) == 0):
        f, ax = plt.subplots(subploty, subplotx, sharex='col', sharey='row')
    # Calculate subplot window index
    ii = (detNum % (subplotx*subploty)) /subplotx
    jj = (detNum % (subplotx*subploty)) % subplotx

    
    for i in xrange(len(dat)):
        if correctGain:
            x = np.arange(256).astype(float) * gainShift[i,detNum]
        else:
            x = np.arange(256).astype(float)
        ax[ii][jj].plot(x, dat[i][:,detNum].astype(float)/dat[2][:,detNum].astype(float) + 1e-10, label = '%s' %filenameList[i])
#        ax[ii][jj].plot(x, dat[i][:,detNum]/acquisitionTime[i] + 1e-10, label = '%s' %filenameList[i])
#        ax[ii][jj].plot(x, dat[i][:,detNum]/acquisitionTime[i]/pulseRate[i] + 1e-10, label = '%s' %filenameList[i])
    
    ax[ii][jj].grid()
    if (ii == 0) & (jj == 0):
        ax[ii][jj].legend()
    #ax[ii][jj].set_xlim((100, 200))
    ax[ii][jj].set_xlim((0, 160))
    
#    ax[ii][ii].set_ylim((1e-2, 1e4))
#    ax[ii][jj].set_yscale('log')
    ax[ii][jj].set_ylabel('Counts')
    ax[ii][jj].set_xlabel('Bin')
    ax[ii][jj].set_title('Det# %d' %(detNum+1))

